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The validation and parallel implementation of a numerical method for the 
r- 1 " 

solution of the time-dependent Dirac equation is presented. This numerical 

method is based on a split operator scheme where the space-time dependence 

is computed in coordinate space using the method of characteristics. Thus, 

o: 

in- 

\o- 

^. ■ spurious solutions related to the fermion-doubling problem and that it can 

be parallelized very efficiently. We consider a few simple physical systems 
such as the time evolution of Gaussian wave packets and the Klein paradox. 
The numerical results obtained are compared to analytical formulas for the 
?-h " validation of the method. 



most of the steps in the splitting are calculated exactly, making for a very 
efficient and unconditionally stable method. We show that it is free from 
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1. Introduction 

The Dirac equation is very important in many fields of physics and chem- 
istry because it gives a relativistic description of electrons and other spin-| 
particles. For this reason, it has been studied and used extensively in rel- 
ativistic heavy ion collisions, heavy ion spectroscopy and more recently, in 
laser-matter interaction (for a review, see [1] and references therein) and 
condensed matter physics [2]. However, solving this equation is still a very 
challenging problem even if it has been derived more than 80 years ago and 
has been utilized profusely. The number of closed-form solutions is very lim- 
ited due to the intricate structure of Dirac matrices which couple the compo- 
nents of the four-spinor wave function. For this reason, only highly symmetric 
systems can be studied by analytical means; the mathematical description 
of more realistic systems should be based on approximation methods such 
as semi-classical theory [3] or numerical calculations. However, the typical 
time scale of the electron dynamics is usually much smaller than the time 
scales of interesting phenomena, rendering the numerical solution notoriously 
difficult and requiring a lot of computer resources. Moreover, certain numer- 
ical techniques (such as the naive symmetric difference scheme) are plagued 
with spurious solutions related to the fermion-doubling problem [4, 5, 6]. Of 
course, such issues should be addressed by any numerical methods to obtain 
accurate solutions describing physically relevant systems. 

The most popular numerical methods rely usually on some variations of 
the split operator method along with a spectral scheme [7, 8, 9, 10, 11]. 
These approaches have also been used to solve the coupled Maxwell-Dirac 
system of equations [12, 13]. In real space, without resorting to a spectral 
method, the finite element scheme [5, 14] and the finite difference scheme 
(both explicit [15] and implicit [16, 17, 18]) have also been exploited. A 



major complication is usually shared by these "real space" approaches: they 
often exhibit spurious states related to the fermion doubling problem [4, 19]. 
The latter leads to the appearance of new unphysical modes when the Dirac 
equation is discretized and is an artifact of the discretization process (a more 
thorough discussion of this subject will be presented in Section 3). To solve 
this problem, some variations of the Wilson [19] and staggered [19, 20, 21, 
22] fermion discretizations have been used. The latter were formulated to 
solve the fermion doubling problem in gauge field theory on the lattice. The 
main issue with these methods is twofold. First, they usually break chiral 
symmetry, one of the fundamental symmetries of the Dirac equation. Second, 
from the numerical point of view, they often lead either to diffusive, unstable 
or non-conservative (probability conservation) schemes. 

In this work, we are implementing and analyzing a variation of the "real 
space" method presented in [23] and show that it circumvents some of the 
issues described previously. This numerical method is a combination of the 
split operator technique and of the method of characteristics; the latter is 
used to obtain 1-D analytical solutions of the Dirac equation in each dimen- 
sion. These solutions are utilized to evolve the wave function at every time 
step by alternate direction iteration and thus, the time-evolution is exact for 
almost every step of the operator splitting (by choosing carefully the time 
and space increments). The resulting scheme is very similar to an upwind 
scheme with a Courant-Friedrichs-Lewy condition CFL = 1. Because exact 
solutions are used in the time evolution, there is no fermion doubling and 
most symmetries are preserved (such as chiral symmetry). Moreover, it yields 
a very simple numerical method which is unconditionally stable and which 
can be efficiently parallelized. 

This article is separated as follows. In Section 2, the numerical method 
and the discretization of the Dirac equation is presented. In Section 3, a 
description of the fermion-doubling problem is exposed along with a proof 



that the numerical scheme is free from these unphysical states. Some details 
of the implementation and code performance, such as the parallelization effi- 
ciency, are shown in Section 4. In Section 5, some simple physical systems are 
analyzed to demonstrate the validity of our numerical approach. The time 
evolution of a Gaussian wave packet and the Klein paradox are evaluated 
and compared to analytical results. We conclude in Section 6. 

2. Numerical Methods 

The Dirac equation describes the relativistic dynamics of spin-| particles 
(fermions) such as electrons and quarks. By construction, it is a one-particle 
theory which is relativistically covariant, i.e. it is invariant under Poincare 
transformations. Other relativistic generalization of the Schrodinger equation 
also exists, such as the Klein-Gordon equation for example. These wave 
functions are distinguished mostly by their spin content (and possibly other 
quantum numbers). For instance, it can be shown, by looking at Lorentz 
transformations, that the Klein-Gordon equation describes spin-0 particles 
while the Dirac equation describes spin-| particles [24]. Of course, the right 
choice of wave equation depends on the system under investigation. 

In this work, we are interested by the relativistic dynamics of an electron 
of mass m coupled to an external classical electromagnetic field characterized 
by its electromagnetic potential. The latter is introduced in the Dirac wave 
equation via the minimal coupling prescription 1 which ensures the invariance 
of the resulting equation under gauge transformations [24]. Therefore, the 
time-dependent Dirac equation is given by [24] 

idti/>(t,x) = Hif>(t,x) (1) 



1 The minimal coupling prescription consists in replacing d^ — > d^ + ieA^, where A^ is 
the electromagnetic potential with Lorentz index [i. 



where ip(t,x) is the time and coordinate dependent four-spinor and H is the 
Hamiltonian operator. The latter is given by 



H = a ■ [cp - eA(£, x)] + (3mc 2 + el 4 V(t, x). 



(2) 



where the momentum operator is p = — iV. Here, A(t, x) represents the 
three space components of the electromagnetic vector potential, V(t,x) = 
Ao(t,x) is the scalar potential, e is the electric charge (obeying e = — |e| for 
an electron), I4 is the 4 by 4 unit matrix and a. = (aij)j,/3 are the Dirac 
matrices. In all calculations, the Dirac representation is used where 



<•■>; 
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Oi oj ' ' ~ [0 -I 2 

The <r,- are the usual 2x2 Pauli matrices defined as 



(3) 
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while I 2 is the 2 by 2 unit matrix. Note that the light velocity c and fermion 
mass m are kept explicit in Eq. (2), allowing to adapt the method easily to 
natural or atomic units (a.u.). 

Throughout this work, we consider the single particle Dirac equation 
which is relevant for calculations describing Quantum Electrodynamics (QED) 
processes coupled to strong classical fields (e.g. particle-antiparticle pair cre- 
ation from very high intensity electromagnetic classical fields). In these con- 
ditions, the negative energy states of the Dirac equation have a well-defined 
physical interpretation and have to be included. However, when QED pro- 
cesses become less important but relativistic effects still need to be taken 
into account, the negative energy states are filled, according to the Dirac 
sea picture. Thus, transitions to the negative energy spectrum are forbid- 
den. Numerically, this physical requirement can be implemented by adding 
projection operators at appropriate places in the calculation (see [17] for a 



discussion of this issue), allowing to consider positive energy states only. This 
however is not investigated in this study, although in principle, it could be 
implemented in our numerical method. 

2.1. Operator Splitting (lowest order) 

The problem at hand is to find an approximate solution of Eqs. (1) and 
(2) at a certain time t n+ i given by ifj n+1 (x) with an initial condition at time 
t n given by 

^(t„,x)=^"(x). (5) 

As previously discussed in [23], this can be done with an operator splitting 
scheme, which forms the basis of our numerical method. To be more specific, 
let us define the operators 



A = 


= — ica x d x 




(6) 


B = 


~ v C-Ct qj (-/•£/ 




(7) 


C = 


= —ica z d z 




(8) 


b = 


= f3mc 2 + eI 4 V(t,x)- 


- ea. ■ A(t,x). 


(9) 



The following splitting in Cartesian coordinates is considered [23] (note that 
we omit the x in the wave functions argument for notational convenience): 

id t ^\t) = Aij (1 \t),^ 1 \t n ) = ij n , te[t n ,t n+1 ) (10) 

id t ^ 2 \t) = B^ 2 \t),^ 2 \t n ) = ^ 1 \t n+1 ),te[t n ,t n+1 ) (ll) 

id t ^\t) = C^(t),^ 3 \t n ) = ^ 2 \t n+1 ),te[t n ,t n+1 ) (12) 

id t ^ A \t) = b^(t),^(t n ) = ^(t n+l ),te[t n ,t n+l ) (13) 

and r +1 = ^\t n+l ) (14) 

where the upper subscript in parenthesis on the wave function denotes the 
splitting step number. Note that this splitting scheme leads to an error that 
scales like (0(5t 2 )), leading to a first order numerical scheme (for more details 
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on the analysis of the method, see [23]). For each step, the initial condition 
and time domain are shown explicitly. Written like this, the method consists 
of solving each equation independently with an initial condition given by the 
solution of the previous step. Note also that for every step, the time incre- 
ment is the same, i.e. St = t n+1 —t n . This way of splitting the Dirac equation 
is very similar to the more usual "T-exponential" operator splitting approach 
used in [7, 8, 13, 12], as can be seen by computing the Fourier transform in 
space coordinates of Eqs. (10) to (12). These numerical techniques, using 
the evolution operator, are usually based on spectral methods which require 
the computation of the discrete Fourier transform at every time step. It is 
possible to circumvent the spectral methods altogether by noting that an 
analytical solutions of Eqs. (10) to (13) can be derived, as discussed in the 
following. 

2. 2. Method of characteristics 

Eqs. (10) to (12) can be solved independently using the method of char- 
acteristics. This essentially proceeds in three steps. First, the Dirac matrix 
is transformed to a new representation where it is diagonalized, thus de- 
coupling the spinor components. The resulting equation has a form similar 
to an advection or transport equation (linear first order in time and space 
derivative). The method of characteristics is then used to find an explicit 
analytical solution. Finally, this solution is transformed back to the original 
Dirac matrix representation. The explicit solution of these equations and 
computational details are relegated to Appendix A. The final result for the 



solutions of Eqs. (10) to (12) is 



/'(Wi-x) = U[h + a x ]^ n (x-c5t,y,z) 
+[I 4 - a x ]i> n (x + cSt,y, z) 



^ (2) (Wi> x ) = -{\^4 + ay\i) {l \t n+ i,x,y - c5t,z) 



(15) 



+[I 4 - ay]ip ( } (t n+1 ,x,y + c8t,z) 



(16) 



1 



* (3) (tn+i,x) = -<{h + a z }ip {2 \t n+1 ,x,y,z-c5t) 
+ [I 4 -a z }ij i2 \t n+1 ,x,y,z + c5t) 



(17) 



Each of these solutions represents traveling waves moving in opposite direc- 
tions at velocity c. 

In the last equation of the splitting, Eq. (13), the solution is simply given 

by 

I'tn+l 



ip { >(t n+l ,x.) = Texp 



x cxp 



rtn+l 

/ dr [(3mc 2 — ea. • A(r, x 



■tn + 1 

—/<: I drV(r : x) 

I'TI 



^\t r 



-1,X) 



:u 



where T, the time-ordering operator, has been introduced. The latter orders 
the argument of the exponential function according to their time argument: 
from the smaller time, on the right to the larger time, on the left. This is 
required because f3mc 2 — ea ■ A(t, x) does not commute with itself when it 
is evaluated at different times due to its Dirac structure (the commutator 
[h(t), h(t')] ^ where h(t) = (5mc 2 — ecu ■ A(t,x)). In the following, this T- 
ordering will be omitted (Texp(...) — > exp (...)) as this approximation results 
in an error of 0(5t 3 ) [7], which has the same (or better) accuracy as the 
method considered. Therefore, this last step of the splitting is the only one 



which is not exact: the time-ordering is neglected and the computation of 
the integration requires a numerical quadrature in general. 

Armed with these analytical solutions, we can now discuss precisely how 
this system of equation is discretized spatially. 

2.3. Spatial Discretization 

In this numerical method, a finite volume formulation is considered where 
the space domain is discretized in cubic elements with edges of length a = 
5x = 5y = Sz. Inside each element, the value of the wave-function is con- 
stant (P type elements). Mathematically, the discretized wave function and 
electromagnetic field can be written as 

TV 

4> h (t,i) = ^i m (i)^(t,x m ) (19) 

m=l 

N 

A h {t,i) = ^l m (i)A(£,x m ) (20) 

m=l 

where N = N x N y N z is the total number of elements, iph{t, i) and Ah(t, i) are 
the discretized wave-function and electromagnetic field (i = (i,j, k) G Z 3 are 
indexing the volumes), the function l m (i) is 1 in volume m indexed by i and 
zero outside, while x m is the vector pointing to the center of volume m. The 
latter is defined as 

1 1 1. \ 

aW + (l + 7)K 3/min + U + 7,)°' ^ min + ( k + 9)° J ( 21 ) 

where x min , y m i n , -^mm are the lower domain boundary coordinates. 

Eqs. (19) and (20) represent the space discretization based on finite 
volume elements of edge size a. On the other hand, Eqs. (15) to (18) are the 
time discretization using an operator splitting with a time increment of St. 
To get the full numerical scheme, they need to be combined together. This 
is performed in the following way. 
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First, note that a priori, the time (St) and space (a) increments are arbi- 
trary and unrelated to each other. However, on the characteristics curves, the 
condition cSt = a is obeyed (this is the equation of the characteristics where 
the partial differential equation becomes an ordinary differential equation). 
This is the reason for the appearance of factors like (x, y, z) ± cSt in the wave 
function argument of Eqs. (15) to (17). Choosing specifically this relation 
between the time and space increments guarantees that the method is un- 
conditionally stable [23]. Moreover, it allows to write ijj(x±cSt) — > iph(i^K) 
when the system of equation is discretized in space and this is exact, i.e. this 
discretized scheme reproduces exactly the analytical form in each dimension 
(up to errors coming from the projection of the grid). This is true for any 
K e N* such that cSt = Ka. In the following, we chose K — 1 for simplicity. 

These facts can be understood in a more physical way. As stated earlier, 
Eqs. (15) to (17) are traveling waves moving at velocity c. Therefore, during 
a time St, the wave travels a distance d = cSt. When the wave is discretized, 
we take its value at the center of each element. By choosing a = d, the value of 
the continuous and discretized wave coincide at every time steps because the 
information moves from one element center to the neighbor element center. 
This would not be the case of other choices of space and time increments and 
this would induce numerical diffusion. 

With this choice of space discretization, Eqs. (15) to (18) become 



t 



m i 



1>. 



n 2 i 



t 



n-.i I 



t 



n+1. 



-{ [I 4 + a x \r h (i ~ 1,3, k) + [h - a x ]M(i + 1,3, k) } (22) 



-\ [I 4 + OyWfrJ ~l,k) + [h - ay]W(i,j + l,k)\ (23) 



-\ [I 4 + a z W(t, k,k-l) + [1 4 - a z W h \i,3, k + 1) } (24) 



cxp 



-i(3mc 2 St - iVh(i) + la. ■ A fc (i) 



t 



n:i I 



(25) 
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where we defined 

pt n +l 

A n (x) = e drA(r,x) (26) 

J t n 

V n (x) = e drV(T,x). (27) 

The discretization of A and V proceeds as in Eq. (19) using the same finite 
volume formulation. The integrals in the last equations can be evaluated 
either analytically, if the expression of the electromagnetic field is simple 
enough, or numerically, using some kind of quadrature. In the latter case, 
this requires a subdivision of the time increment and thus, is more time 
consuming. Moreover, if the Dirac solver is coupled to a Maxwell equation 
solver that computes the time-dependent electromagnetic potential, as in 
[23], the evaluation of these integrals requires even more resources. However, 
if the space-time variations of the electromagnetic field are smaller than the 
typical space-time variations of the wave-function, that is 

|VA| < |W| and \d t A\ < |<9^| (28) 

for all x, then we can simplify Eqs. (26) and (27) as 

A"(x) w A"(x)5t ; V n (x) « V n (x)St. (29) 



This approximation was used in [23] and is valid for many cases of interests 
because the macroscopic electromagnetic field usually involves much larger 
time-scales as that of the electron. In this work however, we consider only 
simple electromagnetic fields which can be integrated analytically and thus, 
the numerical integration is not required. This is also true for certain laser 
pulse models (see [11] for instance for a pulse with a linear turn-on ramp) 
for which the integral can be found analytically. 

The first exponential function in Eq. (25) is the formal expression for a 
4x4 matrix. It could be evaluated numerically by performing a similarity 
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transformation that diagonalizes the matrix in the exponential. However, it is 
more convenient and efficient to use the well-known result for the exponential 
of Dirac matrices, that is 2 

e , (/ 3F +a .F) = co S (| F |) + ^F + a-F) gin(|F|) (3Q) 

\F\ 

where F = (F 1; F%, F 3 ), F, Fi are arbitrary scalar functions and \F\ = \/F 2 + F • F. 
From this result, the solution in Eq. (25) becomes 

r h + \i) = C/(i)exp[-^"(i)]^ 3 (i) (31) 

where C/(i) is a matrix given explicitly by 



U(i) 



•Ah^iX) ( A\ [iAh,x(i)+Ah,y(i) 



C (A) - i^s(A) Mm(!)_M1 s (a) _,^) S (A) 

i^l) s (A) !^i^i s (i) C (A) + i^s(A) 

M h „(i)-A h , y (i)] s(A) _,i^W s(A) c ( A )h-^ s (A) 



(32) 



In the last expression, we defined 



c(A) = cos(A) , s(A) = sin(A) and A = \J (mc 2 5t) 2 + A h (i) • A h (i). (33) 

2.4- Higher order splitting 

The main numerical error in the scheme presented in the preceding section 
is due to the operator splitting. As argued in the last section, this scheme 
is exact in each dimension up to the projection of the wave function and 
electromagnetic field on the mesh and it can be shown that it is a first order 



2 This relation can be derived by using the Taylor expansion of the exponential function 
and the following Dirac matrices properties for any positive integers n: f3 2n — 1, af n = 
1, ^™+ 1 =^a 2 " +1 =a,)- 
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method [23]. This can be improved significantly by choosing higher order 
splitting schemes like the second order Strang-like splitting scheme given by 
(we omit the x in arguments for notational convenience) 



id t ip (6 

id t ip (s 
id t ijj( 9 






B^ 2 

Aifjt 3 

Cy (6 

A^ 7 
A^ 9 



] (t), 


jfjW 


>(*), 


^,(2) 


>(*), 


^,(3) 


>(*), 


ft** 


Ht), 


^,(5) 


>(«), 


^,(6) 


Kt), 


^(7) 


Ht), 


^(8) 


>(*), 


^(9) 



(tn) = 
(tn) = 

(*n+f 
(0 = 
(*„) = 



^ (1) (W±)> ^ 

= ^ (5) (Wi), 
= ^ (6) (Wi), 

= / 8) (Wi), 



£ G 
£ G 
£ G 
£ G 
£ G 
£ G 









(9) 



n+i) m+1, 



n+|) '■n+l, 



and f +1 =^ 9) (t n+1 ) (34) 



Here, the time increments are St/4 = t , 1 — £„ for the operator A, St/2 = 

4 

t t+ i —t n for 5 and C, and finally 5£ = t n+ \ —t n for /}. This kind of splitting 
induces an error of 0(St 3 ). Then, the space discretization proceeds as in 
the lowest order splitting, leading to a second order numerical method [23]. 
Of course, this can be improved to arbitrary order but this increases the 
computation time significantly. In this work, we consider only the lowest 
order and second order splitting. 

It is also important to note that in 1-D and 2-D, the last equation sim- 
plifies considerably because the number of steps in the splitting is reduced. 
For instance, if one considers the x — z plane in 2-D, it is possible to omit 
every steps having the operator B (the same is true for the lowest order 
splitting in Eqs. (10) to (13)) while the subsequent steps with operator A 
can be combined into one step with a time increment of St/2. In 1-D, for the 
x-line, it is possible to omit all steps with operators B and C. In that case, 
the lowest order and second order splitting are equivalent. Also, the Dirac 
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equation in 1-D reduces to a two-component equation and thus, describes the 
time-evolution of a bi-spinor. This simplification however is not implemented 
in the following applications. 

2.5. Time step size 

The time step in the numerical solution of the Dirac equation should be 
chosen carefully to insure the convergence of the calculation. The limit on the 
maximum time step size is related to the fermion mass energy mc 2 because 
the solution (see Eq. (18)) includes a term having a form like ~ e imc t . The 
latter oscillates at a frequency ~ mc 2 , which is very high in comparison to 
other typical time scales. Therefore, the time step should be chosen such as 
[25] 

5t < — - « 5.0 x 10~ 5 a.u.. (35) 

mc 2 

As will be seen in the examples of Section 5, accurate results can be obtained 
with St < 1.0 x 10~ 5 a.u.. 

3. Fermion-Doubling Problem 

To our knowledge, the fermion-doubling problem was first encountered in 
Quantum Field Theory (QFT) on the lattice (see [19, 26] for a review of this 
issue). In these studies, where the goal is the evaluation of path integrals 
for the evaluation of Euclidian correlators, this problem appears when the 
fermionic action (from which the Dirac equation is obtained by a variational 
derivative) is discretized using a symmetric difference scheme for the space 
derivatives 3 . The latter choice is required to keep the hermiticity of the 



3 We define the symmetric difference scheme as d x u(x) —¥ x 2a , where a is the 

lattice spacing. 
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action and of the momentum operator. This however modifies the dispersion 
relation, which becomes (in ID and assuming a free and massless theory) [4] 

E 2 ~ — sin 2 (pa) . (36) 

where a is the lattice spacing, E is the fermion energy and p is the momentum. 
It is then easy to show that the resulting discrete fermionic propagators 
have two poles in the Brillouin zone at p = and p — -, representing two 
different fermionic states (in D dimensions, we have 2 D fermionic states [4]). 
Therefore, the discretized action leads to a theory describing many fermions 
while its continuum counterpart has only one fermionic species. Finding 
a consistent way to resolve this discrepancy between the continuum and 
discrete theories is still an open problem in Lattice QFT. Under fairly general 
assumption, it was proven in the Nielsen-Ninomiya theorem [27] that it is 
actually impossible to keep both the hermiticity of the discretized fermionic 
action and the chiral symmetry without adding the fermion doublers. 

This discrepancy between the continuum and the discretized versions of 
the action also appears at the level of the equation of motion (for coordinate 
space discretizations), as discussed thoroughly in the literature [19, 4, 26]. 
From the QFT lattice studies, two methods have emerged which are free 
from spurious states: the Wilson and the Staggered fermions equation of 
motion. Note however that both breaks chiral symmetry [19]. In the Wilson 
discretization, a new term is added which cancels the second minima of the 
dispersion relation. It is interesting to note that this term is actually the same 
as the artificial viscosity term appearing in advection equation discretization 
schemes (such as in Lax-Friedrichs scheme), which is well-known in numer- 
ical analysis. In the Staggered approach, the fermion doubling is avoided 
by changing the size of the Brillouin zone. This is performed by letting the 
wave function components belong to different lattice points [19]. In general, 
real space discretization schemes of the time-dependent Dirac equation used 
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for applications are variations of these two methods (see [16, 17]), thus cir- 
cumventing the fermion doubling problem. Otherwise, unphysical behavior 
may be seen in the numerical solution. For instance, in a finite element dis- 
cretization of the Dirac equation where the doubling appears, it was shown 
in [5] that an electron wave packet could move at a velocity greater than the 
speed of light. This stresses the importance of this issue. 

We now analyze our numerical method and claim that it does not suffer 
from this fermion-doubling problem. 

3.1. 1-D 

We will now show that our numerical method is free from these spurious 
states. For simplicity, let us first consider the 1-D massless Dirac equation 
without external field (the 3-D case will be treated in the next subsection), 
given by 

id t ip(t,x) = -ica x d x tp(t,x), (37) 

which actually corresponds to the first step of our operator splitting method, 
i.e. Eq. (10). The solution to this equation can be computed by taking the 
Fourier transform in space. From this procedure, we get 

$(t + St, p x ) = e~ ica ^%, Px ) (38) 

where the widehat denotes the Fourier-transformed function and p x is the 
momentum in x. 

Using our discretization method described in the last section, we have 
that 

W +1 (i) = ±{[I 4 + a x ]W(i-l) + [I 4 -a x }W(i + l)}. (39) 

We then take the discrete Fourier transform in space to get 

r h + \k) = l -^[l l + ax }e'^ k + [l A -a x }e^ k )r h {k) 

= e~ ia ^^(k) = e- tca * Pk5t ^(k) (40) 
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where N x is the number of elements, pk is a discrete momentum with index 
k e [0, N x — 1] and the widehat now denotes the discrete Fourier-transformed 
function. To obtain the last equality, we used the fact that in our discretiza- 
tion, we have 5x = c5t. This clearly has the same form as the analytical 
solution Eq. (38) except that the momentum p k = -^^ now takes discrete 
values. Also, we note that by comparing Eqs. (40) and (38), the exponential 
function is just the evolution operator. Therefore, the discrete Hamiltonian 
is Hk = ca x pk from which the energy spectrum can be computed by finding 
its eigenvalues [4]. This yields a dispersion relation given by 

E k = ±cp k (41) 

where E^ is the discrete energy. This is very similar to the continuum case 
except that the momentum and energy take discrete values. The relationship 
between them is linear as in the continuum and therefore, there is no fermion- 
doubling. This result is actually not surprising as our technique is based on 
the exact solution of the Dirac equation. 

3.2. 3-D 

The 3-D case is very similar to the 1-D case. We start by computing 
the analytical solution in the continuum (without mass and without external 
field) by using the Fourier transform method. We get 

$(t + 6t,p) = e - ic{axPx+azPz+azPz)5t i){t^) 



L - iHnSt 



^(t : p) + 0[(5t) 2 ]. (42) 



where Hq is the massless free Hamiltonian (with m = V = A = 0). For the 
discrete part, each coordinate is treated sequentially as in Eq. (10) to (12) 
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because we are using an operator splitting method 4 . Therefore, we get 

^™ +1 (k) = e - ica * Pk > St e- ica « Pk y St e- ica * Pk * 5t fy(k). (43) 

This of course cannot be compared directly to Eq. (42) because of the anti- 
commuting nature of Dirac matrices. However, if we also expand the expo- 
nential, we get an equation similar to the continuum case at the order 0(St), 
allowing us to determine the discrete Hamiltonian and energy spectrum by 
inspection: 

H k = ca z p kz + ca y p ky + ca x p kx + 0(St), (44) 

E^ = ±c^pl+pt y +pl+0(5t). (45) 

This is again a discrete dispersion relation without fermion-doubling problem 
and this proves our assertion that our numerical method is free from these 
spurious states. Note however that in the 3-D case, the dispersion relation 
is not exact as in 1-D due to the operator splitting, which necessitates the 
addition of 0(St) terms and higher. The order of the numerical error on the 
dispersion relation will obviously vary with the splitting order. Neverthe- 
less, this will not change our conclusion concerning fermion-doubling. For 
instance, in the first order splitting, the 0(St) terms in the Hamiltonian are 
H o(St) = c 2 pk ^ ^ a j + c 2 p kx p kz [a z , a x ] + c 2 p ky p kz [a z , a y ] which are linear 
functions of momenta and thus, do not induce fermion doubling. 

4. Implementation and performance 

A computer program using the method described in Section 2 to solve 
the time-dependent Dirac equation was written in C++. Given its relative 
simplicity, it is possible to implement the 1-D, 2-D and 3-D cases easily; 



4 This analysis is for the first order splitting. The analysis can be generalized to higher 
order splitting in the same way. 



changing the number of dimensions requires only minor modifications to the 
code, amounting to the adjustment of some array sizes and to the addition 
(or deletion) of the action of certain operators defined in Eqs. (6) to (8). 
The first order and second order splitting are included; higher order schemes 
could also be easily added but require more computation time and are not 
implemented at the moment. The number of operations per time iteration 
n op scales like n op ~ O(N) where N is the number of elements. This is 
slightly better than spectral methods in which the computation of the Fast 
Fourier Transform scale like n op ~ 0(N In N). In the latter however, there 
are less levels in the splitting scheme because the action of the operator a • p 
is computed in one step instead of three for our numerical method (for lowest 
order splitting). This results in a potentially smaller scaling prefactor. Nev- 
ertheless, the strength of our method is that it can be easily parallelized using 
a domain decomposition strategy, which leads to a very scalable and efficient 
computation tool (see Section 4.1 for the parallel efficiency evaluation). 

4-1. Parallelization 

The algorithm is parallelized by using a domain decomposition strategy: 
the whole domain is divided in subdomains and the mesh data in each of 
these subdomains is sent to a different process. The communication between 
these subdomains, which is required in the computation by the elements 
close to the subdomain boundaries, is performed using the Message Passing 
Interface (MPI). 

To evaluate the efficiency of our parallelization procedure, we look at 
the time-evolution of a 3-D Gaussian wave packet. A mesh of 256 x 256 x 
256 elements is chosen and the wave function is evolved in time for 500 
time iterations (this is the larger problem that can be runned serially using 
our computer resources). The number of processes is varied while keeping 
the computation parameters fixed, thus evaluating the strong scaling. The 
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Figure 1: Parallelization efficiency as a function of the number of processes. 

computation time is recorded for each process number and used to evaluate 
the computation efficiency given by e = -^r- where T p is the computation 
time for p processes. This value, usually between and 1, characterizes 
the quantity of work done by each processes for the actual computation 
in comparison to communication times. For an ideal system with linear 
speedup, we should have e« 1. All the computations are performed on the 
MAMMOUTH cluster 5 . The results for the scaling properties are shown in 
Fig. 1. 

In this figure, it is shown that our numerical method have a super linear 
speedup, i.e. e > 1, for the problem considered (3-D Gaussian wave packet) 



5 The MAMMOUTH is a cluster of Intel Xeon 64 bit, 3.6 GHz CPUs with an Infiniband 
network SDR Cisco- Topsprin non-blocking having a bandwidth of 800 MB/sec. 
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and on the computer used. This at first may seem unexpected: one would 
expect the parallel version efficiency to be smaller than 1 because at each 
time step, processes spend time to communicate with their neighbors, which 
does not happen in serial calculation. The most likely explanation for this 
phenomenon is the cache effect [28]: as the number of processes is increased, 
so is the cache memory (the latter is much faster than ordinary memory). 
Thus, a larger percentage of the problem can be fitted in this memory and this 
decreases the memory access time. The overall effect is to reduce the com- 
putation time significantly. This however implies that the efficiency should 
diminish for larger mesh where this cache effect becomes less important. This 
was observed in calculations using larger mesh where doubling the number 
of processes (for instance, from 64 to 128) resulted approximately in halving 
the computation time (in these cases, a serial computation was not possible). 
Therefore, the efficiency in that case is much closer to 1. 

5. Applications to simple systems 

In this section, some simple systems are studied to validate the numerical 
method and to show its strength. Many simple physical systems are analyzed 
which allow a comparison between analytical formulas and numerical results. 
These are the Gaussian wave packet evolution and the Klein paradox. 

All the computations are performed in atomic units where the electron 
mass is given bym = 1.0 and the speed of light is c = a -1 where a is the fine 
structure constant given by a ~ 1/137.0359895. Note that in these units, we 
also have the electron charge |e| = 1 and the Plank constant H = 1.0. 

In all the examples considered in the following, the wave function at the 
boundaries dQ is set to ip{x., t)\ dn = 0, which corresponds to a Dirichlet 
boundary condition. This of course can induce spurious reflections when the 
wave function reaches the boundary, and this was actually verified empir- 
ically. However, by choosing a large enough domain, this effect is reduced 
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considerably and these reflections can be neglected. We are presently working 
on the implementation of more sophisticated boundary conditions (transpar- 
ent or absorbing). 

5.1. Time Evolution of Gaussian Wave Packet 

Gaussian wave packets in vacuum are among the most simple systems 
that can be studied using a wave equation. For this reason, it has been 
studied extensively in many contexts, whether as a check for the consistency 
of numerical methods [7, 8, 9], to compare with non-relativistic evolution 
[29] or to study the characteristic features of their time-evolution [30, 31, 32]. 
Physically, they represent the localization of a free electron and thus, are very 
important in many applications. In this part of this work, we are analyzing 
the time evolution of massless and massive wave packets in 1-D and 2-D to 
verify and validate our numerical method and its implementation. 

5.1.1. Massless Wave Packet in 2-D 

In this section, the time evolution of a free massless Gaussian wave packet 
is considered. The motivation for studying such a simple systems is twofold. 
First, it is possible to derive analytical solutions (at certain specific space- 
time points) for the time-dependent wave-function in this case, allowing to 
validate and study the numerical method previously developed. Second, 
because it does not have a mass term nor electromagnetic potential, it allows 
us to test the first steps of the splitting and to use a larger grid (the condition 
in Eq. (35) does not have to be fulfilled). 

The starting point of this analysis is the following initial wave packet 
equation in 2-D [24] 



ijj(t = 0,x,z) =N 



e 1 ^) T (46) 
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where M is a normalization constant and A characterizes the Gaussian width. 
This wave packet represents a spin-up massless electron. The analytical 
solution of the free massless Dirac equation with (46) as an initial condition 
can be found for two specific space-time points: 



1. When x = z = 0, we have: 



^(t,0,0) = M 



1 - VSF^e-^erfi (J^ 



where erfi(z) is the imaginary error function 6 . 



(47) 
(48) 



2. When r = \/x 2 + z 2 = ct, we have: 

^(t,r = ct) = AT 2 F 2 Q, 
fo(t,r = ct) = 



1 3 1 1 t 2 

4'2'2' - A2 



ip 3 (t,r = ct) 

ip4,(t,r = ct) 



-Mt 2 sin(0) 2 F 2 



5 7 3 5 t 2 



(49) 
(50) 
(51) 



4'4'2'2 A 2 
-ATt 2 cos(0) 2 F 2 (^;||-^) (52) 



where 2 F 2 is the generalized hypergeometric series 7 and <fi = arctan(^/a;) 
is the polar angle in real space. 

The calculational details necessary to obtain these solutions are relegated to 
Appendix B. 

The results are shown in Fig. 2 for the two points considered and for the 
first wave function component. The domain used in the calculation had a 
size of 40 x 40 a.u., and was discretized in 1000 x 1000 elements using an 



6 The imaginary error function is denned as ern(z) = — j= L e * dt. It is related to 
the error function as erfi(z) = — ier£(iz). 

The generalized hypergeometric series is defined as p F q (ai,...,a p ;bi,...,b q ;x) — 
SfcLo (b)'(b) T\ w ith the Pochhammer symbol defined as (a)k = a(a + l)...(a + k— 1). 
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order 1 splitting. The theoretical and calculated results show a very good 
agreement. 




x \|/ (t,r=0) (calculated) 

i|/ (t,r=0) (theoretical) 

\|/ (t,r=ct) (calculated) 
_\|/ (t,r=ct) (theoretical) 



0.05 0.1 

Time (a.u.) 



0.15 



Figure 2: Time evolution of the first wave function component for a massless wave packet. 

5.1.2. Wave Packet in 1-D 

In this section, the time evolution of a massive wave packet is analyzed. 
The calculation is similar to the one performed in the last section, except that 
the mass term is now considered and the 1-D case is calculated. As discussed 
in Section 2.5, this requires a smaller time step size and consequently, a 
smaller element size (remember that in our numerical method, we have cSt = 
a). The expression of the initial wave-packet is given by Eq. (46) but letting 
z = (ipivit — 0,x) = ip2D(t — 0,x,z — 0)). Again, it represents a spin-up 
massive free electron. 
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An analytical solution for the wave packet time evolution is computed in 
Appendix C and is compared with the result obtained from our numerical 
method. Integrals in Eqs. (C.8) and (C.ll) are calculated in Matlab using a 
very high order integration scheme (adaptive quadrature scheme based on a 
Gauss-Kronrod pair (15th and 7th order formulas)). 

The calculation using our numerical method is performed using the first 
order splitting and for different element/time step size. The domain bound- 
aries are set to ±10 a.u., the wave packet initial position is xq = 0.0 a.u. 
and its width is fixed to A = 1.0 a.u.. We consider three time step size: 
S^ = 3.56 x 10~ 5 a.u., St 2 = 1.78 x 10~ 5 a.u. and St 3 = 8.91 x 10~ 6 a.u., 
corresponding to 4096, 8192 and 16284 elements respectively. The results at 
time t = 2.85 a.u. are shown in Fig. 3. This figure shows how critical the 
choice of time step is, as discussed previously. However, if the condition in 
Eq. (35) is fulfilled, as in the case St 3 , the analytical and numerical results 
are in very good agreement. 

5.1.3. Wave Packet in 2- D 

The same calculation as in the last section is carried in 2-D where the 
initial wave function is given by Eq. (46). The analytical solution is calcu- 
lated in Appendix B and the integrals in Eqs. (B.5), (B.7) and (B.8) are 
again calculated in Matlab using the same numerical method. The domain 
investigated is subdivided into 8192 x 8192 elements and its boundaries are 
positioned at ±5 a.u. (in both x and z). The second order scheme is used, 
making for a time step of St = 1.78 x 10~ 5 a.u.. The wave packet initial po- 
sition is xq = 0.0 and its width is given by A = 1.0 a.u.. The wave function 
is evolved for 32000 time iterations to a final time of t = 0.57 a.u. and the 
results are shown in Fig. 4. The theoretical and calculated results are in 
good agreement for the wave function density p(x, t). This is also true when 
considering the wave function spinorial components as shown in Fig. 5 (for 
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Figure 3: Massive 1-D Gaussian wave packet at t = 2.85 a.u. for different time step sizes. 
The theoretical curve is obtained by integrating numerically Eqs. (C.8) and (C.ll). Note 
that the theoretical result and the third time step (St^ — 8.91 x 10~ 6 ) are overlapping and 
thus, are not distinguishable. 

the first component only). 

5.1.4- Travelling Wave Packet in 1-D 

So far, the calculations performed considered only wave packets at rest. 
In this section, we consider the case where the wave packet is given an initial 
momentum, allowing it to travel on the x-axis. The main reason for studying 
this system is to investigate the phase of the wave function and to determine 
if the latter is reproduced accurately by the numerical method. Initially, the 
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Figure 4: Massive 2-D Gaussian wave packet at t = 0.57 a.u.. Here, the polar coordinate 
r = \/x 2 + z 1 is used because the probability density is symmetric in the azimuthal angle. 
The results shown in the figure corresponds to the positive x-axis (z = and x > 0). The 
calculated results overlap with the theoretical one. 

wavefunction is chosen to be a Gaussian wave packet as 



W 



0,x) 



M 



e ifcox e - 



(x-xq) 

4(A) 2 



with C = 



ck 



me 



+ V 



m 2 c 4 + c 2 k?. 



(53) 



where ko is the wave packet momentum, A is the wave packet spreading and 
xq is its initial position. This corresponds to a superposition of positive and 
negative energy free solutions propagating on the x-axis. 

The analytical solution for the wave packet time evolution is computed in 
Appendix C and is compared with the result obtained from our numerical 
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Figure 5: Massive 2-D Gaussian wave packet at t = 0.57 a.u.. Here, the polar coordinate 
r = \/x 2 + z 2 is used because ipi is symmetric in the azimuthal angle. The results shown 
in the figure corresponds to the positive x-axis (z = and x > 0). The calculated results 
overlap with the theoretical one. 

method. We consider only the first wave function component given by Eq. 
(C.4). The integral in this equation is calculated in Maple using a numerical 
scheme based on an adaptive Gauss-Kronrod quadrature (with Gauss 30- 
point and Kronrod 61-point rules), well suited for oscillatory integrands. 

The calculation using our numerical method is performed using the first 
order splitting and for different element/time step size. The domain bound- 
aries are set to ±5 a.u., the wave packet initial position is xq = 0.0 a.u., its 
width is fixed to A = 0.05 a.u. and its initial momentum to ko = 100 a.u.. 
We consider three time step sizes: 5t\ = 8.91 x 10~ 6 a.u., 5£ 2 — 4.45 x 10 -6 
a.u. and St 3 = 2.23 x 10~ 6 a.u., corresponding to 4096, 8192 and 16284 ele- 
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ments respectively. The results for dh/ji at time t = 0.022 a.u. are shown in 
Fig. 6 (Ssipi and the other components have similar behavior). This figure 
shows that the wave function obtained from the numerical method is not as 
accurate as in the previously studied cases: for the same grid size, the phase 
error is more important. Thus, a smaller time step is required to reproduce 
the theoretical wave function precisely (~4 times smaller). Nevertheless, the 
analytical and numerical results are in very good agreement for St = Sts, 
which confirms the convergence of the numerical scheme. 



Theoretical 

5t = 8.91e-6 

-■- 5t = 4.45e-6 

St = 2.23e-6 




Figure 6: Real part of ipi for a massive 1-D traveling Gaussian wave packet at t = 0.022 
a.u.. The theoretical curve is obtained by integrating numerically Eq. (C.4). Note that 
the theoretical result and the third time step (St^ = 2.23 x 10~ 6 ) are overlapping and thus, 
are not distinguishable. 
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5.2. Klein Paradox 

The "Klein paradox" concerns the reflection and transmission of plane 
wave solution on a potential step of height Vq [33]. In the non-relativistic 
(Schrodinger) case, when the energy of the incident wave is lower than the 
potential step, the solution in region where V(x) is non-zero (in IR + ) is a de- 
caying exponential. As a consequence, most of the wave function is reflected, 
resulting in a very small transmission coefficient. In the relativistic (Dirac) 
case, there is a regime, when E < Vq — mc 2 , where a new phenomenon 
appears: a plane wave solution is now possible, yielding a non-negligible 
transmission coefficient even if E < Vq. This "paradox" has been interpreted 
in the context of the Dirac see picture [34, 35, 36] and in second quantization 
[34, 36] as the production of antimatter on the potential boundary (see also 
[37, 2] for recent views on the subject and possible experimental verifica- 
tion). Thus, the transmitted part is related to the negative energy solution 
which describes anti-fermions, while the incoming and reflected parts are the 
electron wave function. In this section, we consider the scattering on a step 
potential of travelling wave packets in 1-D and 2-D. 

5.2.1. 1-D 

In this subsection, the Klein paradox is analyzed along the same lines 
as in [7, 34], but using the numerical method developed in the preceding 
sections. This numerical test is chosen because it is one of the few existing 
analytical solutions of the Dirac equation for time-dependent systems. The 
calculation is performed in 1-D here and in 2-D in the next section. 

Initially, the wavefunction is chosen to be the traveling Gaussian wave 
packet given in Eq. (53). Also, a smooth potential barrier is considered to 
avoid numerical problems related to discontinuous functions. The potential 
is given by 



V ( 
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x 
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V(x) = — 1 + tanh - (54) 



where Vq is the magnitude of the potential and L characterizes the gradient 
and width of the step. This potential has already been studied for incident 
plane wave solutions and one finds a transmission coefficient given by [34, 38] 

sinh(7rA;L) smh(7ik'L) 
' ~ ~sinh [tt (^ + k + k') §] sinh [n (f - k - k') f] ( ' >j) 

where k = \\/{E k - V ) 2 - m 2 c 4 , k' = -\\/ E 2 - m 2 c A and £ fc = ^k 2 c 2 + m 2 c 4 . 
This transmission coefficient is valid, strictly speaking, for monochromatic 
plane wave solutions. However, the latter are not very convenient for nu- 
merical calculations because they necessitate the addition of sources (and 
absorbers) on the numerical domain boundaries. By choosing A -1 <C fc 
(a Gaussian peaked on momentum k ), the error due to the utilization of a 
Gaussian wave packet instead of a monochromatic wave function is negligible. 

The parameters for the initial wave packet and step potential are the 
same as in [7]: the momentum and width of the wave packet are ko = 106. 
a.u. and A = 1.0 a.u. respectively while the width of the step is L = 10~ 4 . 
The simulation domain is —20 a.u. < x < 20 a.u. which allows us to place 
the initial wave packet at an initial position x = —10 a.u.. For the first 
order splitting, the domain is divided into 65536 elements such that each 
element have a width of a = 6.10 x 10~ 4 (and St = 4.45" 6 a.u.) while for 
second order splitting, we have twice as much elements (131072 elements 
with a = 3.05 x 10~ 4 and the same time increment). 

The height of the potential barrier is varied and the reflection (R) and 
transmission (T) coefficients are measured after a time of 0.22 a.u. which 
corresponds to 50000 time iterations. These coefficients are calculated as 

h+Mx)\ 2 L |^(x)| 2 , x 

The results of the computation for the transmission and reflection coefficients 
are shown in Fig. 8 for the first order splitting while a typical initial and 
reflected wave functions are shown in Fig. 7. 
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Figure 7: Typical wave density (^(x)! 2 ) before and after the interaction of the wave packet 
with the potential barrier. The latter is situated at x = 0.0 a.u.. The transmitted part is 
in the region x > 0, while the initial and reflected wave functions are in x < 0. 

Relative differences between the analytical and numerical results are smaller 
than 0.5% for all values of potential height. Moreover, the transmission co- 
efficient is null for Vq < 4.1 x 10 4 a.u., which is also in accordance with the 
analytical analysis. 

5.2.2. 2-D 

The same system is now analyzed in 2-D. Again, the initial state is a 
Gaussian wave packet prepared at t — given by 



i/j(t = 0,x) =M 



(x-x n ) 2 + (z-z n f 

e lkoX e ^aP with 



(57) 



with a given momentum fc = 100. a.u.. Its width is A = 1.0 a.u. in both x 
and z dimensions. The potential has a height Vq = mc 2 + 1.0 x 10 4 a.u. and 
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Figure 8: Transmission and reflection coefficients for different values of potential step 
height, for the first order splitting (there are similar results for the second order splitting). 
The "calculated" value is obtained from our numerical method while the "theoretical" 
value is obtained from Eq. (55). 
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a thickness of L = 0.0025 a.u.. The domain boundaries are set at x — ±15.0 
a.u. and z = ±5.0 a.u.. The result of the computation for the wave function 
density is shown in Fig. 9 at a time of t — 0.214 a.u., after the wave packet 
has scattered with the potential barrier. In this figure, it is very easy to see 
the transmitted (for x > 0) and reflected parts (x < 0) of the wave function, 
in accordance with the 1-D analysis. The value of the transmission and 
reflection coefficients are 0.2407 and 0.7593 respectively, which is relatively 
close (a relative difference of less than 9% for the transmission and 3% for the 
reflection coefficients) to the analytical results which are 0.2621 and 0.7380 
respectively. 

To obtain this accuracy and to resolve the quasi-discontinuous potential 
barrier, it was necessary to use a very large mesh. The result plotted in Fig. 
9 were obtained with 24576 x 8192 elements and a second order splitting 
(thus, a time step of St = 1.78 x 10~ 5 a.u.). To evolve the wave function to 
the final time t = 0.214 a.u., 12000 time iterations were necessary. With such 
a large mesh, the computation time was approximately T comp . ps 22 hours on 
256 processors 8 , showing the performance of our numerical method. A better 
accuracy would require a mesh with the same element size as in the 1-D case, 
i.e. around a ~ 3 x 10~ 4 (for the same a, it was verified that the 1-D and 2-D 
cases give approximately the same transmission and reflection coefficients). 

6. Conclusion 

In this work, a numerical method for the solution of time- dependent Dirac 
equation in real space was analyzed and used in simple applications. One of 
its main features is that it does not suffer from the fermion-doubling problem, 
as demonstrated in our analysis. The method is based on a combination of the 
operator splitting and characteristic methods and thus, is relatively simple. 



i Again, this calculation was carried on the MAMMOUTH cluster described previously. 
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Figure 9: Final wave function at t = 0.214 a.u. after it scattered with the potential step. 
The transmitted wave function can be found in the region x > while the reflected part 
is in x < 0. Note that the potential step is centred on x = 0. 
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This simplicity resulted in very good parallelization performance and in very 
efficient calculations. Its main drawback however is the relation between 
the time and space increment which has to be chosen as cSt = KSx. This 
condition limits the flexibility for the space grid size in certain circumstances: 
for instance, if the space variations of the wave function are small, it would 
be more efficient to use larger element size but still keeping a small timestep 
(the latter is often necessary because there is a term like e~ tmc * in the Dirac 
equation solution which oscillates rapidly). This is possible with the method 
described in this paper only if some modifications to the numerical scheme 
are implemented. For instance, one could use a Lax-Wendroff scheme to solve 
Eqs. (10) to (12) instead of the method of characteristics, but this would 
introduce again the fermion doubling problem. Note that it is possible to 
analyze systems for which the space variation of the wave function is large 
and the time variation is small by choosing K > 1. 

We applied this method to study simple systems such as the Klein paradox 
in 1-D and 2-D, and the time evolution of Gaussian wave packets (other more 
complex systems are under investigation). The results obtained from our 
method were compared to analytical computations and showed a very good 
agreement. However, to obtain the required accuracy, very small meshes 
had to be utilized, especially for reproducing the phase of the wave function. 
One interesting avenue to circumvent these difficulties would be to consider 
higher order splitting (such as third of fourth order splittings) [39]. Another 
interesting possibility is the use of the GPU architecture which can lead to 
significant improvement in computation time [11]. Work is under way to 
verify if these techniques would result in better computation performance. 
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Appendix A. Solution of the Dirac Equation in Each Dimension 

This Appendix describes how the solution of the Dirac equation in each 
dimension is obtained. We are interested in solving the following equation: 

id t i/)(t,yL) = —icaidiip(t,yi) (A.l) 

with an initial condition given by 

^(t o ,x) = 0(x). (A.2) 

Note here that the index % = x, y, z is not summed and thus, each dimension 
is treated independently. In the first step, the matrix a, is diagonalized to 
decouple the spinor components. This is performed by a similarity transfor- 
mation as 

P}aiPi = A, (A.3) 

where Pi are unitary matrices and Aj = Diagfl, 1, — 1, —1] = /3 is a diagonal 
matrix. Starting with «j in the Dirac representation, the explicit expression 
of the transformation matrix is 

P i = -L(P + a t ) . (A.4) 

The representation thus obtained is very closed to the Majorana representa- 
tion, modulo the index of the matrix in Dirac representation that becomes 

The resulting Dirac equation is then 

id t tp'(t, x) = -icpdiip'(t, x) (A.5) 
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where we defined ■*/>' = P'i/j. It is convenient here, for notational purposes, 
to split the four-spinor into two bi-spinors as ip' = (<//, x') T t° get 

id t (f'(t,x.) = —icdi(p'(t,x.) (A. 6) 

id tX '(t,x) = icd iX '(t,x) (A.7) 

Therefore, the Dirac equation clearly becomes a set of four uncoupled first- 
order differential equations in this representation. Their solutions are well- 
known and can be obtained from the method of characteristics. They are 
given by 

¥>i,2(*> x ) = 4>'i :2 {xi-c5t,x), (A.8) 

xi, 2 (^ x ) = (t>' 3A (xi + cSt,x), (A.9) 

along with the conditions x ± xq = cSt. The latter is the characteristics 
equation along which the partial differential equation becomes an ordinary 
differential equation. Here we also use a particular notation with a function 
argument (£, Xi ± ct,x), which means that the substitution x^ — > Xi ± ct 
is performed while the other coordinates x = (xj,Xk) (with i ^ j ^ k) 
are unchanged. Note that the initial conditions are related to the original 
representation as 0' = P^<p. 

In the last step, the solution is transformed back to the original represen- 
tation. The final result, after some manipulations, is given by 

^(t,x) = -{[I 4 + a i \(f)(xi-c5t,x) + [l4-a i }(f)(x i + c5t,x)}.(A.10) 

This equation is used in the numerical method to evolve the wave function 
in time in alternate direction iteration. 

Appendix B. Solution of 2-D wave packet 

In this Appendix, the analytical solution for the time evolution of 2-D free 
wave packet is computed. We consider an initial wave packet for a massive 
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spin-up electron at rest. The initial wave function is given by Eq. (46) and 
its Fourier transform by 



ilj(t = 0, Px , Pz )=4irA 2 N 



-A 2 (p2+p2) 



(B.l) 



The 2-D Dirac equation we want to solve is given by 

id t tp{t,p x ,p z ) = [ca x p x + ca z p z + [3mc 2 ] ${t,p x ,p z ), (B.2) 

here expressed in Fourier space. The solution to this equation is then simply 



$(t,Px,Pz 



-ica x p x t—ica z p z t—i/3mct 



4>(o,Px,Pz) 

' ftPA .ca x p x + ca z p z + (3mc 2 . ' 

JU cos {hit) — 1 sin {ht) 



E 



xip{0,p x ,p z 



(B.3) 



(B.4) 



where E = \/p x c 2 + p 2 c 2 + m 2 c 4 . This last equation can be Fourier trans- 
formed back to real space. Then, using polar coordinates {p 2 = p 2 + p 2 ) and 
properties of the Bessel functions, we get the solution as 



i)i{t,x,z) = 2AfA 2 / dppe- A ' 2p2 J {pr) 

Jo 



x 




TfiC 

cos(Et) — i—=- sin{Et) 
E 



(B.5) 
ii 2 (t,x,z) = (B.6) 

poo 

i/> 3 (t,x,z) = -2VVA 2 sin(0) / dppj^pr)^ sm{Et)e~ A2p2 (B.7) 

Jo E 

ifj A {t,x,z) = -2A/"A 2 cos(0) / dppJi(pr)^sin(£i)e~ AV (B.8) 

Jo E 

where J n (z) is the Bessel function of the first kind and <fi = arctan(z/:r) is 
the polar angle in real space. There is no known solution for these integrals 
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in the general case [40] but they can be found for two specific cases when the 
mass is zero (m = 0.0): at r = and at r = ct. If r = 0, we have 



^i(t,0,0) -- 


= M 


"i r ct - { -Kt? c( ct \ 
1 — v^r— re ^A^Erfi — - 

v 2A \2Ay_ 


(B.9) 


Mt,o,o) ~- 


= $,(*, o, o) = Mt, o, o) = o 


(B.10) 



where Erfi(z) is the imaginary error function. 
If r = ct we have 



i/>i(t,x,z) 
ip 2 (t,x,z) 
ip 3 (t,x,z) 

ip 4 (t,x,z) 



" % '' - l 4' 4' 2' 2' A 2 







-A/-(ct) 2 sin(0) 2 F 2 t.T^.s; 



5 73 5 
4'4'2'2 : 



-M(ct) 2 cos(0) 2 F 2 



5 7 3 5 
4'4'2'2' 



A 2 

(ct) 2 
A 2 



(B.H) 
(B.12) 
(B.13) 

(B.14) 



where 2 F 2 is the generalized hypergeometric series. These two results will be 
used to validate our numerical method and to analyze the operator splitting. 
These two analytical solutions at the origin (at (x,z) = (0,0)) and at x 2 + 
z 2 = ct for the time evolution of a 2-D massless wave packet are useful mostly 
for the comparison with the numerical method and for validation purposes. 



Appendix C. Solution of 1-D wave packet 

In this Appendix, the analytical solution for the time evolution of a 1-D 
free wave packet is computed. This calculation follows closely the one for 
the 2-D wave packet in Appendix B and for this reason, only the main steps 
are presented. The initial wave function is given by Eq. (53) and its Fourier 
transform by 



$(t = Q,p x ) =2v^r AA/" 



^A 2 ( Px -k ) 2 



(CI) 
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The 1-D Dirac equation we want to solve is given by 

id t $(t,p x ) = [ca x p x + [3mc 2 ] ^(t,p x ), (C.2) 

here expressed in Fourier space. The solution to this equation is then simply 



4>{t,Px 



' ftr ., .ca x p x + /3mc 2 . 

I4 cos (Et) — 1 sin [Et] 

E 



^(0,p x ), (C.3) 



where E = \/p x c 2 + m 2 c 4 . This last equation can be Fourier transformed 
back to real space and we get the solution as 

A 



ip 2 (t,x) 
ip 3 (t,x) 
tp 4 (t,x) 



M 

X 

0, 
0, 

M 

X 



dpe 



-A 2 (p-k ) 2 ipx 



cos(Et) 



, mc 2 ^ cp 



, sm(Et)-iC^sm(Et) 
E E 



A 



dpe' 



-A 2 (p-k ) 2 e ipx 



C cos(Et) + iC— sin(Et) - i— sin(Et) 



(C.4) 

(C.5) 
(C.6) 

(C.7) 



There are no analytical solutions to these integrals, so they have to be com- 
puted numerically. 

We can find the solution for an initial wave packet for a massive spin-up 
electron at rest by setting ko = to obtain 

A 



ip 2 (t,x) 

ip 3 (t,x) 



M 

x 

0, 

0, 



dpe 



-A p ipx 



TflC 

cos(Et) — i—— sm(Et) 

E 



(C.8) 

(C.9) 

(CIO) 



-*A/"A r dpe- A2p2 e ipx —sin(Et). (C.ll) 
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